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MONTE  CARLO  RADIATIVE  TRANSFER  SIMULATIONS  FOR  OCEAN  OPTICS: 

A  PRACTICAL  GUIDE 


1.  INTRODUCTION 

This  report  provides  a  methodology  and  the  equations  needed  to  generate  Monte  Carlo  computer 
simulations  for  many  common  ocean  optics  applications.  These  applications  include  the  modeling  of 
natural  ocean-atmosphere  environments  and  analyses  of  laboratory  and  in-situ  optical  instrumentation. 
We  attempt  to  provide  enough  practical  detail  to  make  it  straightforward  for  the  reader  to  write  his/her 
own  computer  code  for  his/her  own  application.  This  document  also  serves  as  documentation  of  the 
methods  the  authors  have  used  in  several  specific  research  projects  (e.g.,  [1] — [3]). 

Monte  Carlo  simulations  in  this  context  involve  the  tracing  of  the  fates  of  millions  of  virtual  photons 
or  photon  packets  according  to  statistical  probabilities.  The  source  of  the  photons  being  simulated  can  be 
either  the  sun  or  an  artificial  light.  Each  simulated  photon  path  is  randomly  distinct  from  the  others  as 
determined  by  the  probabilities  of  absorption  and  scattering  in  the  water  and  air  and  by  the  probabilities 
of  transmission,  reflection,  or  absorption  at  air-liquid  and  liquid-solid  interfaces.  Example  ray-tracing 
diagrams  of  Monte  Carlo  simulations  are  shown  in  Figure  1  for  the  open  ocean  water  column  and  for  a 
point-source  integrating  cavity  absorption  meter  (PSICAM).  In  the  ocean  case,  photons  enter  at  the  sea 
surface,  scatter  within  the  water,  reflect  internally  at  the  sea  surface  or  seafloor,  and  eventually  either 
escape  into  the  atmosphere  or  are  absorbed  by  the  water  or  seafloor.  In  the  PSICAM  case,  photons  are 
initiated  at  the  center  of  the  cavity,  internally  scattered  by  the  water,  internally  reflected  by  the  cavity 
wall,  and  eventually  absorbed  by  the  water  or  cavity  wall.  The  light  intensity  at  any  region  in  these 
modeled  systems  is  determined  by  counting  the  number  of  photon  paths  that  intersect  that  region. 


Figure  1  -  Forward  Monte  Carlo  simulations  of  an  ocean  water  column  (left)  and  of  an  integrating-cavity  absorption  meter 

(right). 
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The  Monte  Carlo  approach  provides  the  most  general  and  most  flexible  technique  for  numerically 
solving  the  equations  of  radiative  transfer.  Because  Monte  Carlo  techniques  are  computationally 
expensive,  they  have  historically  been  avoided  whenever  alternative  solution  techniques  are  available. 
However,  as  computer  speeds  continue  to  dramatically  improve,  the  Monte  Carlo  approach  is  becoming 
more  practical  and  popular  for  all  types  of  radiative  transfer  calculations.  Furthermore,  it  is  no  longer 
imperative  that  Monte  Carlo  codes  be  complicated  with  tricks  designed  to  save  computer  memory  and  to 
avoid  expensive  commands  (i.e.,  logical  statements  and  trigonometric  functions).  This  makes  Monte 
Carlo  techniques  accessible  now  to  anyone  with  basic  programming  skills. 

Monte  Carlo  techniques  have  been  used  by  various  researchers  in  the  ocean  optics  community  over 
the  last  several  decades  [1] — [20],  and  introductions  to  the  use  of  Monte  Carlo  techniques  in  ocean  optics 
are  provided  in  Refs.  21  and  22.  References  on  the  use  of  Monte  Carlo  methods  for  the  study  of  light 
propagation  in  the  atmosphere  [23]  and  in  biological  tissue  [24,  25]  are  also  relevant.  Most  of  the  cited 
references  discuss  very  specific  applications  and  fail  to  provide  sufficient  detail  to  reproduce  their  Monte 
Carlo  simulations;  however,  they  at  least  provide  insight  into  general  Monte  Carlo  methodology.  In  this 
report  we  compile  much  of  the  practical  information  provided  in  the  references  and  add  lessons  learned 
from  our  own  research.  We  introduce  some  basic  ocean  optics  and  Monte  Carlo  concepts  in  Section  2, 
describe  how  to  move  a  photon  from  one  interaction  point  to  another  in  Section  3,  and  show  how  to 
compute  scattering  angles  in  Sections  4.  A  description  of  how  to  simulate  photon  sources  and  detectors  is 
provided  in  Sections  5  and  6.  The  interactions  of  photons  with  surfaces,  such  as  the  sea  surface  and 
seafloor  are  discussed  in  Section  7.  Backward  Monte  Carlo  techniques  are  introduced  in  Section  8.  A 
short  summary  of  the  most  frequently  used  Monte  Carlo  equations  is  provided  in  Section  9.  In  Section  10 
we  provide  many  examples  of  Monte  Carlo  codes  relevant  to  Ocean  Optics  and  evaluate  their 
performance. 


2.  PRELIMINARIES 
2.1  Nomenclature 

The  terminology  in  this  report  is  a  combination  of  that  from  radiative  transfer  theory  (as  applied  in 
optical  oceanography)  and  that  from  probability  and  statistics.  For  the  former,  we  adopt  the  nomenclature 
from  Mobley  [21].  For  example,  the  processes  of  absorption  and  scattering  by  seawater  are  quantified  by 
the  beam  absorption  coefficient  a  (m1),  the  scattering  coefficient  b  (m1),  and  the  scattering  phase 
function  [:i  .  The  beam  attenuation  coefficient  c  (m1)  and  single  scattering  albedo  ay  are  related  to  a  and  b 

by  c  =  a+b  and  ay  =  b!c.  These  symbol  definitions,  as  well  as  others  used  in  this  report,  are  provided  in 
Table  1. 


Table  1  -  Nomenclature 


A 

beam  absoiption  coefficient  (m'1) 

B 

beam  scattering  coefficient  (nf 1 ) 

& ri  br;q 

Raman  scattering  coefficient 

c 

beam  attenuation  coefficient  (m1),  c  =  a+b 

E0 

scalar  irradiance  (W/m2) 

Ed 

downward  irradiance  (W/m2) 

Eu 

upward  irradiance  (W/m2) 
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G 

scattering  phase  function  asymmetry  factor 

i,j,k 

unit  vectors  along  x,  y,  and  z  axes 

L 

optical  pathlength 

N 

index  of  refraction 

n 

outward  normal  to  a  surface  (unit  vector) 

P 

cumulative  distribution  function  (CDF) 

P 

probability  density  function  (PDF) 

R 

irradiance  ratio 

Rb 

seafloor  albedo 

Ri 

internal  reflectance  at  water-air  interface 

q 

random  number  uniform  on  [0, 1  ] 

s 

geometric  pathlength  (m) 

w 

photon  weight 

x,  y,  z 

Cartesian  coordinates 

X 

position  vector,  x  =  [xi ,  yj,  zk \ 

Zb 

water  depth  (m) 

p 

scattering  phase  function 

e 

polar  angle  with  respect  to  the  z-axis 

9 

incident  polar  angle  with  respect  to  the  normal  to  a  surface 

9, 

transmitted  polar  angle  with  respect  to  the  normal  to  a  surface 

X 

Wavelength 

V 

Frequency 

M* 

cosine  of  polar  angle  with  respect  to  photon  direction  or  with 
respect  to  the  normal  to  a  surface. 

Mx 

direction  cosine  with  respect  to  x-axis 

My 

direction  cosine  with  respect  to  y-axis 

Mz 

direction  cosine  with  respect  to  z-axis 

p 

surface  reflectivity 

T 

optical  depth 

azimuthal  angle  with  respect  to  direction  of  photon  travel 

</> 

azimuthal  angle  with  respect  to  the  x-axis 

1 

unit-magnitude  direction  vector 

y 

polar  angle  with  respect  to  direction  of  photon  travel 

Q 

solid  angle 

a>0 

single-scattering  albedo,  a>Q=blc 

2.2  Probability  Functions  and  Random  Numbers 

Monte  Carlo  simulations  rely  on  properly  defined  probability  density  functions  and  cumulative 
distribution  functions.  In  general  the  probability  density  function  p(x )  is  defined  such  that  the  probability 
that  a  specified  event  will  occur  over  the  range  x  to  (x+cLr)  is  \p(x)  dx]  and 

poo 

p(x)dx  =  1. 


(2.1) 
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The  cumulative  distribution  function  P(x)  is  the  probability  that  the  specified  event  will  occur  between  the 
lowest  possible  value  of  x  (frequently  0,  but  more  generally  -oo)  and  x, 


P(x)  =  f  p(x)dx,  0  <  P(x)  <  1 . 

J  —oo 


(2.2) 


To  determine  the  value  of  x  for  a  particular  Monte  Carlo  realization,  we  select  a  random  number  q 
from  the  uniform  distribution  over  [0,1],  set  it  equal  to  the  cumulative  distribution  function  [21], 

P(x)  =  q,  (2.3) 

and  solve  for  x.  We  solve  for  x  either  analytically  or  numerically,  depending  on  the  form  of  P(x).  Many 
examples  of  the  use  of  Eq.  (2.3)  are  provided  throughout  this  report.  It  is  important  to  note  that  a 
particular  value  of  q  is  never  used  more  than  once.  If  consecutive  equations  contain  q,  each  q  represents  a 
distinct  random  number.  This  is  consistent  with  the  way  random  number  generators  are  used  in  computer 
languages.  For  example,  when  writing  code  for  Matlab  (The  MathWorks,  Inc),  the  variable  q  can  simply 
be  replaced  with  the  keyword  rand.  At  each  and  every  occurrence  of  rand,  Matlab  will  substitute  a  new 
random  number. 


2.3  Photon  Weights  and  Biasing 


Because  it  is  computationally  wasteful  to  trace  the  full  paths  of  photons  that  never  reach  a  spatial 
region  of  interest,  it  is  frequently  prudent  to  sample  from  a  biased  cumulative  distribution  function. 
Biasing  the  distribution  function  makes  it  possible  to  trace  more  photons  that  are  likely  to  find  their  way 
to  areas  of  interest  and  fewer  photons  that  are  unlikely  to,  without  changing  the  final  computed  result. 
This  makes  it  possible  to  converge  more  quickly  to  a  result  of  specified  accuracy.  If  the  probability 
function  is  p(x ),  but  we  prefer  to  trace  photons  according  to  the  biased  probability  density  function  pi,(x) 
so  that  more  photons  reach  regions  of  interest,  we  can  use  pdx)  as  long  as  we  multiply  the  photon  weight 
by  the  ratio  of  the  correct  distribution  to  the  biased  one  [21], 


w-w 


'  P(x)  " 

yPb(X), 


(2.4) 


An  example  of  biased  sampling  is  given  in  Sect  3.1;  additional  examples  are  provided  by  Mobley  [21]. 

We  use  a  similar  idea  to  treat  the  absoiption  of  photons  within  the  water  column  and  at  surfaces. 
Rather  than  computing  the  probability  of  absoiption  and  then  terminating  the  photon  path  if  the  photon  is 
absorbed,  we  trace  a  photon  to  an  event  (scattering  or  reflecting)  and  then  reduce  its  weight  based  by  the 
statistical  probability  that  absorption  would  have  taken  place.  We  then  continue  to  trace  the  photon  path 
as  if  it  had  not  been  absorbed.  An  alternative  way  to  think  of  this  is  that  we  initially  launch  a  packet  of 
many  photons  traveling  in  the  same  direction  and  we  reduce  the  number  of  photons  in  the  packet  at  each 
event  by  the  number  of  photons  that  would  have  been  absorbed  by  the  event.  Specifically,  we  multiply  the 
photon  weight  by  the  value  of  the  single-scattering  albedo  a>o  each  time  it  is  scattered  and  by  the  value  of 
the  surface  albedo  each  time  it  is  reflected  by  a  surface.  Only  when  the  weight  of  a  photon  is  reduced  to 
below  a  certain  specified  level  do  we  stop  tracing  the  photon  path. 
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3.  PHOTON  TRAVEL 
3.1  Path  Length 

From  the  definition  of  optical  distance  /,  the  probability  density  function  for  the  attenuation  of  light 
with  respect  to  optical  distance  traveled  is  [21] 

p(l)  =  e~l,  l>  0.  (3.1) 

From  Eq.  (2.2),  the  cumulative  distribution  function  is  [21] 

/ 

P(l)  =  \e-rdl'  =  \-el  .  (3.2) 

0 

To  determine  /  for  Monte  Carlo  simulations,  let  P(f)=  q  [Eq.  (2.3)]  and  solve  for  /  [21], 

/  =  -ln(l- q)  =  -ln(^) ,  0<^<1.  (3.3) 

Shown  in  Figure  2  are  plots  of  p(l),  P(l),  and  l(q)  as  given  by  Eqs.  (3.1)  -  (3.3). 
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0  0.2  0.4  0.6  0.8  1 


q 


Figure  2  -  Probability  functions  for  pathlength  determination  and  pathlength  (l)  versus  random  number  (q). 

In  homogenous  waters,  the  geometric  pathlength  5  (in  meters)  can  be  computed  with  [21,  22] 

_  l  _  ln(g) 
c  c 


(3.4) 


where  c  is  the  attenuation  coefficient  (m'1).  In  nonhomogenous  waters,  the  relationship  between  s  and  /  is 
given  by 


/  =  ]"  c(s')ds' ,  (3.5) 

and  the  value  of  s  must  be  computed  by  starting  at  zero  and  increasing  its  value  until  the  total  value  of  / 
computed  with  Eq.  (3.5)  equals  that  given  by  Eq.  (3.3).  In  a  layered  system,  the  integral  in  Eq.  (3.5) 
becomes  a  summation. 


A  simple  Monte  Carlo  code  that  relies  primarily  on  Eq.  (3.4)  is  provided  in  Section  10.1. 
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Equation  (3.3)  gives  the  probability  density  function  for  the  distance  an  unimpeded  photon  travels 
before  being  scattered  or  absorbed.  However,  if  there  is  a  highly  absorbing  boundary  present,  it  is 
computationally  wasteful  to  track  a  high  number  photons  that  collide  with  this  boundary.  We  can  instead 
sample  from  a  biased  density  function  that  only  considers  photons  that  are  scattered  before  reaching  the 
boundary  [19], 


_eip(-0_ 

l-exp(-/6) 


o  <  /  <  4 , 


where  lh  is  the  optical  distance  to  the  boundary  along  the  direction  of  travel.  To  compensate  for  the  biased 
probability  function,  we  multiply  the  weight  of  the  photon  packet  by  [1  -  exp(-/b)],  which  is  the  fraction  of 
the  photons  that  do  not  reach  the  boundary.  The  cumulative  distribution  function  is 


P(l)  =  j‘0P(Odl' 


Q  l  +  lb  _  g  l  1 

-elb  +  2-  e~'b 


Note  that  P(lh)  =  1.  Solving  P(l)  =  q,  we  obtain 

/  =  4-ln[exp(4)  +  (l-exp(4))^]. 


3.2  Direction  Cosines 

The  direction  cosines  //x,  jUy,  and  //z  are  the  x,  y,  and  z  components  of  the  unit  vector  that  points  in  the 
direction  of  photon  travel.  By  definition,  they  satisfy 

ju;  +  ju;  +  /4  =  I . 


Given  the  direction  cosines  and  the  photon  path  length  s,  we  can  advance  the  photon  from  its  previous 
position  (xo,yo,zo)  to  its  next  location  with 


X  =  x0  +  jUx  s  , 


y  =  y  o  +  My  s  > 

Z  =  Z0+/JzS. 

For  a  direction  defined  by  polar  angle  0(0  <6  <7t)  and  azimuthal  angle  tj)  ( 0  <  tj)  <  2n ), 

0 <0<tt,  0<(j)<2n , 


the  corresponding  direction  cosines  are 


fix  ^sin^cos^, 

(3.6) 

/j  =sin^sin^, 

(3.7) 

jLL  =  cos  6  . 

(3.8) 

Leathers,  Downes,  Davis,  and  Mobley 


Because  0  <  6<  n,  sin^is  always  positive  and  can  be  determined  from  //z  with 

sin6>  =  ^l  -  ful  . 

Frequently  we  know  //z  and  (f>,  from  which  we  can  find  //x  and  //y  with 

M,  =  V1-/'--  cos^’ 

My=yl1~/J=  sin^- 


To  compute  6  and  <j>  from  the  direction  cosines  we  can  use 

0  =  cos-1  juz , 


>  =  cos 


/4 


K 


Mv  ^ 0 


>  =  2^- -cos  1 


A 


■//; 


,  //_:  *  l,  //.,  <  o . 


If  ///  =  1,  then  jux  =  jUy  =  0,  making  <j>  indeterminate.  In  this  case,  <j>  can  be  set  randomly. 

3.3  Time-Dependent  Simulations 

If  we  wish  to  compute  how  long  it  takes  for  light  to  find  its  way  from  source  to  detector, 
compute  the  time  it  takes  each  photon  to  trace  its  path.  The  time  required  for  a  photon  to 
geometric  pathlength  5  is 


t  =  s/c„  -  sn! cn 


we  can 
move  a 


where  Co  is  the  speed  of  light  in  vacuo  and  cn  is  the  speed  of  light  in  a  medium  with  index  of  refraction  n. 
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4  SCATTERING 

4.1  Scattering  Probabilities 

4.1.1  Basic  equations 

The  scattering  phase  function /((T,®)  is  the  probability  density  function  that  gives  the  probability 
that  a  photon,  when  scattered,  will  scatter  at  polar  angle  ¥  and  azimuthal  angle  O  away  from  the  incident 
direction  [21],  From  Eq.  (2.1),  the  integral  of  /)( T.®)  over  all  directions  is  unity, 

j  Jo  PQV,  O )  sinCT^ W®  =  1 .  (4.1) 

The  sin(T)  term  in  this  equation  arises  from  the  definition  of  the  spherical  coordinate  system.  For 
seawater  and  for  air,  the  scattered  azimuthal  angle  ®  with  respect  to  the  incident  direction  is  uniformly 
distributed  over  [0,2rc].  Therefore  p(x V)  and  /;(®)  are  independent  of  one  another  and 

pOF,®)  =p(¥)p(&). 


In  order  to  satisfy  Eq.  (2.1), 


r  2  k 

j0  /’(®)o,®  =  l, 


the  probability  density  function  for  ®  must  be 


^(O)  = 


2  n 


(4.2) 


From  Eq.  (2.2) 


fO  1  ® 


In 


and  from  Eq.  (2.3)  [21] 


®  =  2/r<7.  (4.3) 

Alternatively,  one  can  use  two  random  numbers  q\  and  qi_  to  determine  ®  with  [23] 

Wx=\-  2qx ,  W2  =1-  2q2 , 

d  -  sjwf^wf , 


cos  ®  -  Wf  d  ,  sin  ®  =  W2  /  d  . 


(4.4) 
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Because  for  seawater  and  air  /?('F,cp)  does  not  depend  on  ®,  Eq.  (4.1)  simplifies  to 

2^J'  /?(¥)  sin  (¥)*/¥  =  1 .  (4.5) 

From  Eqs.  (2.1)  and  (4.5)  [21], 

pi'V)  =  2;zj0(vP)sin(vP)  .  (4.6) 

From  Eqs.  (2.2),  (2.3),  and  (4.6)  [21], 

P(x F)  =  p(x¥)dx¥  =  fiC¥) sin('F)fi?vP  =  q.  (4.7) 

Alternatively,  if  we  express  the  scattering  phase  function  in  terms  of  the  cosine  of  the  scattering  angle, 

jus  =  cos'F,  dps  =  -sin('F)c/'F 

then  Eqs.  (4.5)  -  (4.7)  become 


-2/rJj  fi(p,)dMs  =  2/rj  j , 6(ps)dps  =  1 , 

(4.8) 

P(M,)  =  2nfi(p,), 

(4.9) 

P(//s)=[  p(ju)dju  =  2x\  fi(ju)dju  =  q. 

"Ms  "Ms 

(4.10) 

For  a  given  phase  function,  one  must  evaluate  the  integral  in  Eq.  (4.7)  or  (4.10)  and  solve  for  T  or  //s  in 
terms  of  q. 

We  note  that  some  authors  define  P(ps)  such  that  its  integral  equals  unity,  as  opposed  to  Eqs.  (4.5) 
and  (4.8)  that  contain  a  2rc  term.  If  so,  then  Eqs.  (4.9)  and  (4. 10)  must  also  have  the  2 n  term  removed;  i.e., 

I"  ~P(p)dp  =  q 

J  Ms 

if  and  only  if 

M)dp  =  \. 

4.1.2  Isotropic  scattering 

For  all  specific  cases  we  will  consider  the  azimuthal  scattering  angle  ®  is  determined  with  Eq.  (4.3), 

®  =  2  nq  . 

For  isotropic  scattering  fi  is  constant  by  definition.  In  order  to  satisfy  Eqs.  (4.5)  and  (4.8) , 
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4  n 

From  Eqs.  (4.6)  and  (4.9), 

f>(T)  =  lsin(T), 


and  from  Eq.  (4.10), 


|  J*'  si n(vE)(/vE  =  |(1  -  cos'?)  =  ? , 
or 

Thus  the  polar  scattering  angle  is  determined  with 

Ms  =  cos  *FS  =\  —  2q. 


Plots  of  /Js(q)  and  'V(q)  for  isotropic  scattering  are  shown  in  Figure  4. 

4.1.3  Rayleigh ,  Raman,  and  pure  water  scattering 

Rayleigh,  Raman,  and  Pure-water  scattering  phase  function  are  of  the  form 

/?ccl  +  />2, 

where / is  a  scalar  constant.  Specifically,  scattering  for  pure  water  obeys  [21,22] 

/?  oc  1  +  0.835//2 , 


(4.11) 


and  Raman  scattering  at  490  nm  follows  [21,22] 

(3  cc  1  +  0.55//2 . 


For  a  given  value  of f  (e.g.,  0.835  for  pure  water)  we  can  normalize  the  probability  density  function 
such  that  its  integral  from  p  =  -1  to  1  is  unity, 


p(m)  = 


3  i+jy 

2  3  +  / 
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Integrating,  we  obtain  the  cumulative  distribution  function, 


P(J*)  = 


f 


6  +  2/ 


M 


3  1 

- /U  H - . 

6  +  2/  2 


To  find  a  value  for  /j  in  Monte  Carlo  simulations,  we  set  P(ju)  =  q  and  solve  for  q.  Therefore,  the  value  of 
/j  is  the  real  root  of 


/ 

6  +  2/ 


M 


3 

6  +  2/ 


u-\ - q  -  0  . 

2 


(4.12) 


Equation  (4.12)  can  be  solved  either  analytically  (which  has  a  complicated  form)  or  numerically.  Figure  3 
shows  the  relationship  between  random  number  q  and  polar  scattering  direction  /us  for  pure  water  as 
computed  by  substituting/=  0.835  into  Eq.  (4.12). 


0.9 

0.8 

. 

0.7 

0.6 

0.5 

0.4 

0.3 

0.2 

0.1 

0 

-1 

-0.5 

0 

0.5 

1 

<? 

Figure  3  -  Determination  of  the  cosine  of  the  scattering  angle  from  random  number  q  for  pure  water. 


4.1.4  Henyey-Greenstein  scattering  phase  function 

The  Henyey-Greenstein  phase  function  [26]  is  defined  in  terms  of  the  asymmetry  factor  g, 

1  1-g2 


An  (1  +  g2  -  2g  cos  ¥  s ) 


3/2 


-Kg<l, 


Pig',  Ps) 


1  1-g 

4^  (l  +  g2-2g/03'2’ 


-Kg<l, 


(4.13) 


where  g  =  0  reduces  to  isotropic  scattering  and  g  close  to  1  describes  very  forward  scattering.  Equation 
(4.13)  satisfies  Eq.  (4.8),  and  from  Eq.  (4.9), 


(l-g2)sin¥, 

(1  +  g2  ~  2g  cos  T/  )3/2  ’ 
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p(ms)  = 


1  1-g2 

2  (1  +  g2  ~  2g/us)312 


From  Eq.  (4.10), 


P(M,)  = 


1-g2 


r 

j  u 


Ms  (1  +  g2  -  2g//.s)3/2 


1-g2 

f 

1 

\ 

1 

2g 

l1_g  1 

Jl  +  g2  -  2gMs  , 

■q- 


Solving  for  /js , 


2g  +  l-2qg-2q  +  2q2g  +  g2  -Ig’q-lqg2  +2giq2 
(~g-\  +  2qg)2 


For  convenience,  this  can  be  rewritten  as 


Ms 


2g 


f 


1  +  g2- 


1-g2 


K\+g-2gRJ 

Ms  =  1  -  2  g, 


g  *0, 

g  =  o. 


Note  that  if  we  substitute  ( 1  -q)  for  q,  we  obtain  an  alternative  form  [24,25], 


Ms 


2g 


1  +  g2- 


1-g2 


l-g  +  2g^ 


g^O  . 


The  value  of  /js  is  plotted  versus  random  number  q  for  several  different  values  of  g  in  Figure  4.  For  g  =  0, 
the  scattering  is  isotropic  (Sect.4.1.2). 
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0  0.2  0.4  0.6  0.8  1 

q 

Figure  4  -  Polar  scattering  angle  versus  random  number  q  for  isotropic  scattering  (g  =  0)  and  Henyey-Greenstein  scattering  with 

anisotropy  factors  of  g  =  0.5  and  0.9 

4.1.5  Measured  scattering  phase  functions 

Few  measurements  of  the  scattering  phase  function  have  been  made  in  natural  waters.  The  most  well 
known  are  those  made  by  Petzold  [27];  other  early  measurements  are  summarized  by  Jerlov  [28].  One 
way  to  use  such  measurements  in  Monte  Carlo  codes  is  to  tabulate  the  cumulative  distribution  function  as 
a  function  of  q  and  use  interpolation  to  compute  T'  for  each  q. 

Kirk  [4]  computed  the  cumulative  distribution  function  for  the  Petzold  measurement  made  in  San 
Diego  Flarbor.  This  is  included  here  in  Table  2  and  Figure  5.  Similarly,  the  Petzold  average  particle  phase 
function  [11]  can  be  integrated  to  obtain  its  CDF. 


Table  2-  Cumulative  distribution  function  versus  scattering  angles  for  Petzold’s  San  Diego  scattering  phase  function  (computed 

from  Kirk  [9]). 


P(0s) 

a 

Ms 

P{0s) 

Os 

Ms 

0.517 

2.5 

0.999048 

0.9832 

92.5 

-0.04362 

0.666 

7.5 

0.991445 

0.9853 

97.5 

-0.13053 

0.751 

12.5 

0.976296 

0.987 

102.5 

-0.21644 

0.811 

17.5 

0.953717 

0.9888 

107.5 

-0.30071 

0.851 

22.5 

0.92388 

0.9902 

112.5 

-0.38268 

0.879 

27.5 

0.887011 

0.9917 

117.5 

-0.46175 

0.901 

32.5 

0.843391 

0.9928 

122.5 

-0.5373 

0.918 

37.5 

0.793353 

0.9941 

127.5 

-0.60876 

0.931 

42.5 

0.737277 

0.995 

132.5 

-0.67559 

0.942 

47.5 

0.67559 

0.996 

137.5 

-0.73728 

0.95 

52.5 

0.608761 

0.9968 

142.5 

-0.79335 
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0.957 

57.5 

0.5373 

0.963 

62.5 

0.461749 

0.968 

67.5 

0.382683 

0.972 

72.5 

0.300706 

0.975 

77.5 

0.21644 

0.978 

82.5 

0.130526 

0.981 

87.5 

0.043619 

0.9976 

147.5 

-0.84339 

0.9982 

152.5 

-0.88701 

0.9988 

157.5 

-0.92388 

0.9992 

162.5 

-0.95372 

0.9997 

167.5 

-0.9763 

0.9998 

172.5 

-0.99144 

1.0000 

177.5 

-0.99905 
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Figure  5  -  Cosine  of  polar  scattering  angle  versus  random  number  q  for  Petzold’s  San  Diego  Harbor  scattering  phase  function. 

To  avoid  having  to  interpolate  tables,  Haltrin  [8]  suggested  using  an  empirical  fits  to  these  functions. 
He  provides  mathematical  expressions  for  many  different  measured  scattering  phase  functions,  including 
those  obtained  by  Petzold,  Kopelevich,  Man’kovsky,  and  others. 

4.2  Updating  The  Direction  Cosines 

Given  an  initial  direction  defined  by  0  and  (f>  and  the  scattering  angles  T  and  ®  with  respect  to  the 
initial  direction,  we  need  to  determine  the  new  direction  cosines,  juf,  juf ,  and  ///.  If  we  let  a  represent  the 
unit  vector  in  the  direction  of  the  initial  photon  direction, 

a  =  [Mx’Mx’Mx] , 


then  the  new  direction  unit  vector  is  given  by 

a  =  sin  T*  cos  ®a±  x  a  +  sin  T*  sin  ®a±  +  cos  'Fa  , 

where  a±  is  a  unit  vector  perpendicular  to  a.  Because  a±  is  not  unique,  there  is  not  a  unique  formula  for 
updating  the  direction  cosines.  If  we  choose  a  to  lie  in  the  x-y  plane,  then 

a±=±[-My,Mx,0]l^~M2z  , 
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a±  x  a  =  ± 


MXM: ,  MyMz  >-(!-//_-)]/  V1-# 


If  the  photon  is  very  close  to  the  z-axis  (e.g.,  |//z|  >  0.99999)  then  it  is  preferable  to  update  the  direction 
cosines  with  [25] 


jUx  ’  -  sin  Y  cos  ®  , 


//,, 1  =  sinxFsin®, 


M-_ 


-!-E-  cos  'F  —  sign(//_ )  cos  T. 

I  I 


Putting  this  all  together,  the  new  direction  cosines  can  be  determined  from  the  initial  photon  direction,  the 
cosine  of  the  polar  scattering  angle  //s  (//s  =  cosT),  and  the  azimuthal  scattering  angle  ®  with  [1,  23,  25], 


/4 

A 

A 


MxM:/yj 

1  -Ml 

-Mv/^-mI 

Mx 

V 

- T 

1  -  ms  cos® 

MyM-Jsl 

'l  -Ml 

Mx/^-mI 

My 

M 

k- A  sin® 

-yfT- 

Ml 

0 

M. 

Ms 

M'x 

- mI  cos® 

My 

ii 

x/l  ~ Ml  sin® 

A  _ 

Ms 

A  <i 


(4.14) 


In  plane  parallel  problems  where  //_-  is  the  only  direction  cosine  that  we  need  to  keep  track  of,  we  can 
update  ik  simply  with  [4] 


/4 '  =  M,_  Ms  +  V1  “ Ml  V1  “ m]  cos®  . 
As  a  consistency  check,  any  transformation  should  satisfy 


Ms  =  YMx,My,MA\Mx\My\Mzr\ 


and 


J{Mxf+{Myf+(M:?=l- 


4.3  Inelastic  Scattering 

In  inelastic  scattering,  some  of  the  light  at  a  given  wavelength  is  removed  and  reemitted  at  a  longer 
wavelength.  Specifically,  we  are  concerned  with  the  Raman  scattering  of  liquid  water.  The  difference  in 
the  frequencies  of  the  incident  and  scattered  light  is  defined  as  the  Raman  shift  vr  [21,  29,  30], 
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v,.=vs  -v\ 

The  amount  of  Raman  scattering  present  can  be  quantified  by  the  Raman  scattering  coefficient  br:q  (nf1), 
which  is  the  ratio  of  scattered  quanta  per  second  at  vs  to  the  incident  quanta  per  second  at  v'.  This 
coefficient  is  an  inherent  optical  property  and  [30] 

a+b+ b  = c , 

although  usually  the  Raman  coefficient  is  implicitly  included  in  the  value  of  a  [21]. 

The  quantity  br  q  is  related  to  br,  which  is  the  ratio  of  scattered  power  at  vv  to  the  incident  power  at  v',  by 
[30] 


K  = 


V' 


\V.J 


The  wavelength  dependence  of  br  is  given  by  [2 1  ] 


bMY 


K 

VT'y 


br{K)> 


where  A0  is  a  reference  wavelength;  br  can  be  taken  to  be  2.4  x  10"4  m"1  for  Xf  =  488  nm. 

One  way  to  handle  inelastic  scattering  in  Monte  Carlo  codes  is  to  simulate  detectors  counting  photons 
at  multiple  wavelengths  (the  excitation  wavelength  and  several  other  longer  wavelengths).  Following 
Waters  [29],  we  used  a  random  number  to  determine  whether  or  not  a  photon  had  undergone  an  inelastic 
scatter  (based  on  the  probability  quantified  by  b,).  If  so,  another  random  number  was  used  to  determine 
the  frequency  shift.  If  the  new  wavelength  for  the  photon  was  within  one  of  the  wavelength  bands  being 
measured,  then  we  continue  to  trace  the  photon  at  the  new  wavelength. 

The  Raman  scattering  direction  can  be  determined  by  substituting /=  0.55  into  Eq.  (4.12), 

-  0.1549///  -  0.4225/4  +  0.5  =  q . 


One  must  either  solve  for  jus  for  each  random  number  q  or  tabulate  //  versus  q  and  then  read  p  from  the 
table.  Example  values  are  provided  in  Table  3. 


Table  3  -  Cosine  of  Raman  scattering  direction  versus  random  number  q. 


q 

Ms 

0.0000 

-1.0000 

0.1111 

-0.8195 

0.2222 

-0.6148 

0.3333 

-0.3841 

0.4444 

-0.1311 

0.5556 

0.1311 
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0.6667 

0.3841 

0.7778 

0.6148 

0.8889 

0.8195 

1.0000 

1.0000 

The  shift  in  wavenumber  can  be  approximated  with  [29]:  a  zero  probability  for  shifts  less  than  2975 
cm'1,  a  linear  increase  from  2975  to  3175  cm'1,  a  uniform  distribution  from  3175  to  3425  cm1,  a  linear 
decrease  from  3425  to  3725  cm'1,  and  a  zero  probability  above  3725  cm'1.  Normalizing  this  piecewise- 
linear  probability  density  function  (PDF),  we  obtain  the  four-point  function  shown  in  Table  4.  A  more 
complex  model  for  the  Raman  scattering  PDF  is  provided  by  Mobley  [21]. 


Table  4  -  Probability  density  function  for  Raman  shift  in  wavenumber  [20], 


wavenumber 
k  (cm1) 

/>(*) 

2975 

0 

3175 

0.002 

3425 

0.002 

3725 

0 

The  cumulative  distribution  function  (CDF)  is  obtained  by  integrating  the  probability  density 
function.  For  practical  use  it  can  either  be  tabulated  or  approximated  by  an  analytic  function.  The  CDF  for 
the  PDF  of  Table  4  is  plotted  in  Figure  6  and  tabulated  in  Table  5. 


K 


Figure  6  -  Cumulative  distribution  function  for  the  Raman  wavenumber  shift  as  computed  from  Table  3. 
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Table  5  -  Cumulative  distribution  function  for  Raman  wavenumber  shift  (from  Figure  6). 


q  =  P(k) 

K 

q  =  P(K) 

K 

0.0000 

2975 

0.6000 

3375 

0.0125 

3025 

0.7000 

3425 

0.0500 

3075 

0.7917 

3475 

0.1125 

3125 

0.8667 

3525 

0.2000 

3175 

0.9250 

3575 

0.3000 

3225 

0.9667 

3625 

0.4000 

3275 

0.9917 

3675 

0.5000 

3325 

1.0000 

3725 

Given  the  wavenumber  shift  A  a;  the  shift  in  wavelength  for  excitation  wavelength  /,  is  [21] 


AA  = 


a 


Ate' 

To7, 


where  a:  is  in  ent  and  A  is  in  nm. 

Shown  in  Figure  7  are  the  probability  functions  of  emitted  Raman  radiation  versus  wavelength  for 
incident  radiation  at  four  particular  wavelengths.  Shown  are  both  the  results  for  the  approximate  PDF  of 
Table  1  and  for  the  function  provided  in  Ref.  21.  It  can  be  seen  that  the  emitted  Raman  radiation  is  spread 
over  a  larger  range  of  wavelengths  for  high  values  of  incident  wavelength  than  it  is  for  small  values  of 
incident  wavelength. 
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Figure  7  -  Raman  emission  probability  functions  for  400,  450,  500,  and  550  nm.  The  piecewise-linear  profiles  are  computed 
from  Table  4,  and  the  smooth  functions  are  reproduced  from  Mobley  [21]. 


5.  SOURCES 


5.1  Determining  The  Initial  Photon  Direction 

5.1.1  Source  Coordinate  System 


To  determine  the  initial  direction  cosines  for  a  photon,  we  first  determine  the  polar  and  azimuthal 
angles  0S  and  </>s  of  the  photon  relative  to  the  frame  of  reference  of  the  source  (or  to  the  normal  to  the 
surface)  and  then  convert  0,  and  (j>s  to  //„  juy,  and  //,.  For  the  special  case  of  a  source  described  in  the  same 
frame  of  reference  as  the  geometry  of  the  problem  (e.g.,  diffuse  illumination  at  a  horizontal  sea  surface), 
then  9=  9S  and  (f>  =  <f>s,  and  the  direction  cosines  are  computed  with  Eqs.  (3.6)  -  (3.8).  If  the  source  frame 
of  reference  is  not  aligned  with  the  axes  of  the  Earth,  then  one  must  compute  the  direction  of  the  source’s 
z-axis  with  respect  to  the  Earth’s  frame  of  reference.  If  we  align  the  source  v-axis  along  the  Earth’s  x-y 
plane,  then  the  initial  direction  cosines  of  the  photon  are  related  to  the  direction  cosines  of  the  source  z- 
axis  /4„  jLiys,  and  jU:s  by  [From  Eq.  (4.14)] 


Mx 

MXSM-_S/ 

'SSt 

~Vys/ 

'SSL, 

Mxs 

sin  0  cos  O 

My 

= 

MysMzs  / 

'SSI, 

Mxs/ 

sl^Ml 

Mys 

sin  0  sin  O 

_/E_ 

-V1 

-  mI 

0 

Mzs 

COS0 

Mx 

sin©  cos  O 

My 

-  sign  ( M-s ) 

sin  0  cos  O 

_Mz. 

COS0 

Monte  Carlo  Radiative  Transfer  for  Ocean  Optics 


21 


where  ©  and  ®  define  the  direction  of  the  photon  with  respect  to  the  z-axis  of  the  source. 

Let  7(0,®)  represent  the  source  function,  normalized  such  that 

j  |  7(0,®)  sin0i/@i/®  =  1 . 

For  azimuthally  independent  sources, 

2k  j  7(0,®)  sin 0<70  =  2^|  I(jus)dju  =  1 .  (5.1) 

For  azimuthally  independent  sources  (or  for  modeling  azimuthally  independent  or  azimuthally 
integrated  quantities)  the  probability  distribution  function  for  the  azimuthal  angle  is  [analogous  to  Eq. 
(4.3)], 

®  =  2  nq.  (5.2) 

The  probability  density  function  for  the  polar  angle  is  therefore 

p(&)  =  2k  7(0)  sin0,  0 <  0 < tt/ 2 , 
or 

7 KM, )  =  )  >  0  ^  Ms  £  1 .  where  ps  =  cos  ©  .  (5.3) 

To  determine  the  initial  photon  direction  in  Monte  Carlo  simulations  we  set  the  cumulative  distribution 
function  equal  to  a  random  number  q, 

P{&)  =  2k j  ° 7(0)  sin0i/0  =  g,  O<0<^-/2, 

P(MS  )  =  27rf  1  (M,  )d/i,=q,  0<Ms<\,  (5.4) 

J/'o 

and  solve  for  0  or  //s. 

For  internal  sources,  let  1(9,  f)  represent  the  internal  source  function  normalized  such  that 

f  [  l(9,(jf)  sin 6  dOd(j>  =  f  [  I(9,</>)  d pdtj)  =  1 . 

J0»0  J  0  J  —  1 

As  with  surface  sources,  the  azimuthal  angle  is  determined  with 

(/>  -2k  q . 

For  azimuthally  independent  problems, 

2k  j  1(9, <j>)  sin  9  d9-  2k  j  I(ju )  dju  =  l. 


(5.5) 
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which  allows  us  to  define  the  probability  density  function  as 

p(6)  =  In  1(0)  sinO, 
pip)  =  2n  I(p) . 

To  determine  the  initial  photon  direction  in  Monte  Carlo  simulations  we  set  the  cumulative  distribution 
function  equal  to  a  random  number  q, 

P(6)  =  2njo °  1(6)  s 'mO  d6-q,  0<6<n!2 
P(p)-2n\  I(p)dp  =  q,  -\<p<\ 

J  rt) 

and  solve  for  0  or  //s. 

5.1.2  Ideal  Collimated  Light 

For  ideal  collimated  light  (e.g.,  a  laser  beam  or  direct  sunlight)  we  can  simply  initiate  all  photons  at 
the  same  location  and  pointed  in  the  same  direction,  all  with  a  weight  of  unity.  For  an  actual  source, 
however,  some  amount  of  variability  should  be  included  in  the  source  location  and  direction. 

5.1.3  Isotropic  Point  Source 

For  an  isotropic  point  source  /  is  constant.  To  satisfy  Eq.  (5.5), 

/  (ds)  =  I(ps)  =  ^-, 

4  n 

so  [22] 

/?(<?,)  =  -  sin<9s,  p(ps)=X-. 

The  cumulative  distribution  function  is 

P(0S )  =  r,  {,/'  sin  Od 0,  P(ps )  =  ^  J'  dp. 

So,  analogous  to  the  result  for  isotropic  scattering, 

cos  0s  =  ps  =1-2  q. 

5.4.3  Diffuse  Illumination 

For  diffuse  illumination,  the  radiance  is  constant  at  all  angles.  To  satisfy  Eq.  (5.1), 
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From  Eq.  (5.4), 


Thus, 


M/O  =  ~  Ms  ■ 
n 


f  1  2 

P(M,)  =  2 1  jud]U  =  \-Ms  =q. 

Ms  =  > 


or  simply 


Ms  =  ■ 


5.2  Initial  Photon  Position 

5.2.1  Ideal  Sources 

For  ideal  sources,  one  simply  initiates  each  photon  at  the  center  of  the  source  location. 

5.2.2  Gaussian  Laser  Beam 

For  a  Gaussian  laser  beam  for  which  b  is  the  He  radius,  the  profile  of  relative  irradiance  E(r)  [mm'2] 
as  a  function  of  radial  position  r  is 


£(,)  = 


exp(-r2/&2) 
nb 2 


The  probability  density  function  describing  the  beam  profile  as  a  function  of  radial  position  r  is  [23] 


P(r)  = 


exp  (-r2/b2) 


2  r,  where 


J  p(r)dr  = 


=  1. 


The  cumulative  distribution  function  P(r)  is 


P(r )  =  | p^/'jdr'  =  I  - 


exp 


f  r^ 
\  b'  5 


Setting  P(r)  equal  to  a  random  number  q  and  solving  for  r  yields 


r  =  b  y]-ln(\-q)  , 
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which  is  then  the  rule  for  sampling  p(r). 

5.2.3  Sources  for  One-Dimensional  Problems 

Many  ocean-atmosphere  problems  can  be  reduced  to  one-dimension,  as  the  desired  results  frequently 
depend  only  on  ocean  depth.  In  this  case,  the  3-dimensional  source  must  be  collapsed  down  to  one 
horizontal  location  and  expressed  in  terms  of  intensity  per  square  meter. 

5.5  Biased  Sampling  Of  The  Source  Function 

5. 5. 1  Upward  biasing  of  source  polar  angle 

Suppose  we  have  an  isotropic  point  source  and  a  detector  located  directly  above  the  source  [21],  As 
shown  above,  for  an  isotropic  source  p(0)  =  sin  612.  Because  we  would  like  to  trace  more  photons  in  the 
upward  direction  than  in  the  downward  direction,  we  can  instead,  for  example,  use 


Pb(0)  = 


\l\  —  £2 

7l(\  +  £  COS  O') 


and  assign  initial  photon  weights  of 


7i(l  +  s cos 6)  sin# 


5.5.2  Biasing  the  source  azimuthal  angle 

Biasing  of  the  azimuthal  angle  must  be  treated  differently  than  polar  angle  biasing.  The  azimuthal 
angle  is  defined  over  [0,2tc],  whereas  the  polar  angle  is  defined  over  [0,Tt].  Also,  the  probability 
distribution  function  for  the  polar  angle  contains  a  sine  term  arising  from  spherical  geometry,  whereas 
that  of  the  azimuthal  angle  does  not. 

If  the  source  is  azimuthally  symmetric,  then  the  proper  probability  distribution  function  is 

=  b)<(j)<2n. 

2  n 

Suppose  we  wish  to  bias  the  photons  toward  <j>=  0,  which  is  along  the  x-axis.  One  function  that 
accomplishes  this  over  [0,2n]  is 


Ph(</>)  = 


\J  \-£2 

2n{\-£  cos^)’ 


D<(j)<2n . 


However,  because  inverse  trigonemetric  functions  return  values  on  [0,7i],  it  is  more  convenient  to  work 
with 


P  >,(</>)  = 


7t(\-£  COS(/>) 


0 <(j)<n , 


Monte  Carlo  Radiative  Transfer  for  Ocean  Optics 


25 


and  make  adjustments  at  the  end.  The  proper  probability  distribution  function  over  0  >  </>>  te  is  (1/ji),  so 
we  assign  the  photons  the  initial  weight 


(1  -  £  COS  0 ) 
sk~£2 


The  corresponding  cumulative  distribution  function  is 


P(#)  =  2 


arctan 


(l  +  f)tan(^2) 

>/C l-ffjO  +  ff) 


/r 


V(i-ff)(i+ff) 


A 


y 


The  resulting  relationship  between  (/)  and  q  is 


=  2  arctan 


tan 


qxyjjl -£)(!  +  £) 

2s!\-£2 


^l(\-£)(\  +  £) 


\  +  S 


J 


If  the  problem  itself  (i.e.,  not  just  the  source)  is  symmetrical  about  the  ^=0  direction  (i.e.,  the  x- 
axis),  we  can  stop  here  and  let  the  photons  all  be  initiated  on  0  >  </>>  n.  Otherwise,  generate  an  additional 
random  number  and  let 


</>  =  <!>,  (q~  1/2)  >0, 


(j> -2k -(f),  — 1/2)  <  0. 
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6.  DETECTORS 

6.1  Detectors  for  3-dimensional  problems 

Generally  we  want  to  measure  either  the  light  that  is  absorbed  by  the  water  or  the  light  field  that 
remains.  The  former  is  easier  to  implement  in  Monte  Carlo  codes;  each  time  a  photon  is  absorbed  we  note 
its  location  and  increment  the  counter  for  that  location  by  the  absorbed  photon's  weight.  To  measure  the 
light  field,  on  the  other  hand,  we  need  to  use  logical  statements  to  determine  if  each  photon  path  crosses 
the  area  of  the  sensor  in  an  appropriate  direction. 

In  three-dimensional  problems  we  generally  treat  Monte  Carlo  photons  as  power  [W],  The  downward 
and  upward  irradiance  measurements  are  computed  with 

4  4) 00  ~lN^Wi’ 

where  A  is  the  area  of  the  detector  and  N  is  the  number  of  photons  traced.  For  radiance  detectors  we 
collect  the  photon  weights  that  reach  the  detector  within  the  appropriate  solid  angle  and  divide  the  result 
by  the  sensor  area  and  solid  angle  Q.  By  definition,  d£2  =  sint?  dt?d^.  For  a  sensor  with  a  conical  field  of 
view  (FOV)  with  half-angle  Of,  the  solid  angle  of  its  FOV  is 

ef 

Q  =  2n  |  sin  Odd  -  2/r(l  -  cos  6f ) . 

0 


For  example,  the  upwelling  radiance  is 


44) 


i 

2/r(l  -  juf)  NA  ’ 


40,  ^  Mf , 


where  ///is  the  cosine  of  the  half-angle  FOV  (///  =  cos6f).  The  radiance  and  irradiance  at  non- vertical 
orientations  are  computed  in  the  same  manner  except  that  the  logical  check  for  admitting  photons  changes 
to  a  comparison  between  juz  and  the  FOV  about  the  normal  to  the  detector. 

6.2  Detectors  for  One-Dimensional  Problems 

The  ocean  is  frequently  approximated  as  an  infinitely  wide,  plane-parallel  medium  in  which  the  only 
variation  is  with  depth.  In  this  case  we  perform  Monte  Carlo  simulations  purely  along  the  z-axis.  The 
Monte  Carlo  photons  in  this  one-dimensional  arrangement  represent  energy  per  area  [W/m2]  rather  than 
energy  [W],  and  we  can  take  the  irradiance  [W/m2]  at  a  given  depth  to  be  proportional  to  the  weighted 
number  of  Monte  Carlo  photons  that  cross  the  xy  plane  at  that  depth.  Taking  positive  z  to  be  the  vertically 
downward  direction  into  the  ocean,  the  downward  and  upward  irradiances  normalized  to  the  input  flux  at 
the  surface  are  given  by 
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4iO)  =  ^2>,.,  (/O,*  °> 

ou<0> 

where  V  is  the  number  of  photons  traced  and  w,  and  (//z),  are  the  weights  and  z-direction  cosines  of  the 
photons  that  cross  depth  z.  Measurement  of  the  irradiance  at  directions  not  perpendicular  to  the  z-axis 
require  that  the  sum  of  the  photon  weights  be  divided  by  the  z-direction  cosine  of  the  normal  to  the 
detector  to  account  for  the  change  in  photon  concentration  (photons  per  area  perpendicular  to  the 
direction  of  travel), 


m = 


i 


Ncosdi 


w, . 


For  an  ideal  scalar  (spherical)  irradiance  detector,  the  detector  surface  is  always  perpendicular  to  the 
travel  direction  of  photons  striking  it,  and  therefore  each  photon  weight  should  be  divided  by  its  own 
direction  cosine  For  example,  the  downward  scalar  irradiance  for  the  one-dimensional  geometry  is 
computed  with 


Radiance  measurements  (W/m2/sr)  are  made  by  counting  only  photons  that  fall  within  the  specified 
field  of  view  (FOV)  and  dividing  the  final  result  by  the  corresponding  solid  angle.  For  the  upwelling 
radiance  the  sensor  is  aligned  with  the  z  axis  and 


40) 


1 

2^(1  -jUf) 


Hi  ^  H  f , 


where  ///  is  the  cosine  of  the  FOV  half-angle.  More  generally,  radiance  at  a  nominal  viewing  angle  0  is 
given  by  [4] 


m= 


Ed(0) 


In  sin0  cos#  Ad 


where  Ed(0)  is  the  irradiance  reaching  the  plane  of  the  sensor  within  the  sensor  FOV. 
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7  SURFACE  INTERACTIONS 

7.1  Air- Water  Interface 

7.1.1  On  the  Air  Side 

When  photons  strike  an  air-water  interface,  a  fraction  of  them  will  be  reflected  and  the  rest  will  be 
transmitted.  We  therefore  need  to  compute  the  fraction  that  is  reflected  and  the  directions  of  the  reflected 
and  transmitted  packets.  The  angle  of  incidence  with  respect  to  the  normal  to  the  surface  h  is 

0i  -  cos-1  \jUz-h\. 

The  reflected  angle  6,  (with  respect  to  the  normal  to  the  surface)  is  the  same  as  the  incident  angle, 

e,=ei. 


and  the  transmitted  angle  6,  is  [21] 


9,  =sin 


—sin  0 


where  n,  and  n,  are  the  indices  of  refraction  for  the  incident  side  (water)  and  for  the  transmitted  side  (air) 
of  the  interface. 

For  the  special  case  of  a  horizontal  surface, 


9i  =  cos-1 


r 

N 

•  -1 

n,  [ 

1  ..2 

sin 

— V 

\-n: 

U 

) 

The  fractional  reflectance  for  unpolarized  light  is  given  by  [21,  22] 


sin(t 

%-<>,) 

2 

+ 

tan(t 

%-o 

) 

1 

sin(< 

+ 

tan  ( ( 

%  +  o 

)J 

1 

p(e,)-- 


Kn{+n2j 


,  e,  =  o. 


If  we  wish  to  model  polarization  effects,  we  can  use 


ni  cos  9i  -  nt  cos  ()t 
nj  cos  Gi  +  n,  cos  Gt 


pMA) 


(7.1) 
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pMA) 


/  \  2 

nt  cos  9t  -  iii  cos  9t 

,  n,  cos  9 i  +  n,  cos  9t  J 


(7.2) 


p{e,A) 


Pv+P\\ 

2 


(7.3) 


To  trace  both  photons  that  are  reflected  and  those  that  are  transmitted,  we  must  determine  each  time  a 
photon  reaches  the  surface  whether  it  is  reflected  or  transmitted.  We  can  do  this  simply  by  drawing  a 
random  number  q  and  letting  the  photon  reflect  if  and  only  if  q  <  p{  9„  9,)  and  transmit  if  and  only  if  q  > 
pi  ty,  Of  If,  on  the  other  hand,  we  are  not  interested  in  the  photons  that  are  reflected  off  the  water,  we  can 
treat  all  photons  as  if  they  are  transmitted  and  multiply  the  photon  weight  by  [1  -  p(9„9t)\.  In  either  case, 
the  total  optical  pathlength  traveled  by  the  photon  from  one  scattering  point  to  the  next  should  equal  the 
value  of  /  computed  with  Eq.  (3.3). 

7.1.2  On  the  Water  Side 

The  angle  of  incidence  with  respect  to  the  normal  to  the  surface  ft  is 

9j  =  cos-1  \juz-h\. 


The  reflected  angle  0,  (with  respect  to  the  normal  to  the  surface)  is  the  same  as  the  incident  angle, 


Or=Oi, 


and  the  transmitted  angle  0,  is  [21] 


f . 


9,  -  sin  1 


—sin  <9 


) 


where  nt  and  n,  are  the  indices  of  refraction  for  the  incident  side  (water)  and  for  the  transmitted  side  (air) 
of  the  interface.  For  the  special  case  of  a  horizontal  seasurface, 

9,  =  cos-1  juz. 


9t  -  sin  1 


\n: 


_ h 


Because  n,  >  n,  (i.e.,  traveling  from  water  to  air),  there  is  a  critical  incident  angle  0C  above  which 
there  is  1 00%  reflection, 


9  =sin 


\ntJ 


The  fraction  reflected  for  unpolarized  light  is  given  by  [2 1  ] 
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p(e,A)=\ 


sin  (0,-0,) 
sin((9;.  +Ot) 


+ 


tan (Oi  -9,) 
tan((7  +Ot) 


Oi<Oc,Qi*  0, 


p{th>n2)  = 


^  -  n2  ^ 


V«i  +  nij 


,  0t  =  O 


P(0„9t)  =  1,  0;>0C. 

Reflectance  for  polarized  light  is  provided  in  Eqs.  (7.1)  -  (7.3).  If  the  surface  is  horizontal,  we  can  check 
if  the  photon  is  beyond  the  critical  angle  by  comparing  juz  to  jit, =cos(  9r), 


If  we  wish  to  trace  both  photons  that  are  reflected  and  those  that  are  transmitted,  then  we  must 
determine  each  time  a  photon  reaches  the  surface  whether  or  not  it  is  internally  reflected.  If  it  arrives  at  an 
angle  greater  than  the  critical  angle  then  it  reflects.  If  not,  we  can  simply  draw  a  random  number  q  and  let 
the  photon  internally  reflect  if  and  only  if  q  <  fi  d,/),).  If,  on  the  other  hand,  we  are  not  interested  in  the 
photons  that  are  transmitted,  we  can  treat  all  photons  as  if  they  are  internally  reflected  and  multiply  the 
photon  weight  by  p{6„9,).  In  either  case,  the  total  optical  pathlength  traveled  by  the  photon  from  one 
scattering  point  to  the  next  should  equal  the  value  of  /  computed  with  Eq.  (3.3). 

7.1.3  Wind-Blown  Surfaces 

Let  9„  and  (j)n  represent  the  polar  and  azimuthal  angles  of  the  normal  to  a  wave  facet  with  respect  to 
the  normal  to  the  horizontal  surface  (i.e.,  Ogives  the  wave  slope,  where  0„  =  0  is  a  horizontal  surface,  and 
(j)n  gives  the  wave  orientation). 

The  cumulative  distribution  function  fort),  is  [20] 


p{0*)  =  l- exp 


— '-r  tan 2  6 
2cj2 


(7.4) 


a2  =0.003  +  0.00512(7. 


From  (7.4)  and  (2.3), 


tan  0n  =  cr^/— 2  In  (l  -q), 

or,  since  q  is  evenly  distributed  on  the  interval  [0,1],  we  can  write 

tan  9n  =  cTy]-2\n{q). 


The  value  of  <A,  is  drawn  from  a  different  random  number  [20,1 1], 
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<t>n  =  2  nq- 


The  slopes  of  the  wave  facet  along  the  x  and  y  directions  are 

dz 

z*=~r  =  tan^  cos^. 

ax 


and 


z 


y 


—  =  tan<9„sin  (j>n. 
ay 


Since  the  unit  vector  normal  to  the  wave  facet  can  be  written  as 

[Z*’V] 


77  = 


2  ,  2  ,  , 

Z,-  +  Z„  +1 


the  angle  of  the  incident  photon  relative  to  the  face  of  the  wave  facet  is 

COS  6i  -  [//,. ,  /Jy ,  |/4  |]  •  77. 


Proper  implementation  of  these  equations  requires  some  additional  logic.  At  the  very  least,  one  must 
verify  that  the  photon  is  not  impinging  on  the  wave  facet  from  the  wrong  side  [20].  Accommodations 
should  also  be  made  for  the  possibility  of  a  photon  interacting  with  more  than  one  wave  in  the  same 
surface  interaction  and  for  the  fact  that  not  all  wave  slopes  are  possible  for  a  given  photon  incident 
direction  [11,21], 

7.2  Seafloor 

7.2.1  General  Equations 

The  interaction  of  light  with  the  seafloor  is  quantified  by  the  bi-directional  reflectance  function  that  is 
a  function  of  both  incident  and  reflected  directions.  Conceptually,  we  like  to  think  about  a  light  beam 
traveling  in  a  particular  direction  (8„  (/>,)  being  reflected  into  another  particular  direction  ( 0r,  </>,).  But 
since  any  source  has  some  finite  divergence,  and  any  detector  has  some  finite  field  of  view,  we  can 
associate  small  solid  angles  dQx  and  df2r  with  the  incident  and  reflected  beams,  respectively.  Let  Lx{0„ 
</>,),  and  L,(  0r,  (/>,.)  denote  the  radiance  of  the  incident  beam  and  the  reflected  radiance  (Figure  8). 
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Figure  8  -  Quantities  used  in  the  definition  of  the  BRDF. 

Our  goal  is  to  define  a  quantity  that  tells  us  how  the  reflective  properties  of  the  surface  vary  with  incident 
and  reflected  directions.  Therefore,  consider  a  measurement  in  which  we  hold  the  direction  of  the 
detector  in  Figure  8  constant  while  we  vary  the  direction  of  the  source.  The  BRDF  is  then  defined  as 


BRDF{6t,4t,erjr) 


dLr{OrA) 

Li(0i,</>i)cos0i  d£lt{pt,fc) 


[sr1]. 


(7.5) 


Note  that  if  we  change  only  the  magnitude  of  the  incident  radiance,  the  reflected  radiance  will  change 
proportionately,  and  the  BRDF  will  remain  unchanged.  However,  if  we  change  the  direction  of  the 
incident  or  reflected  beams  while  holding  all  else  constant,  the  BRDF  will  change. 

Suppose  we  want  to  compute  the  total  radiance  heading  upward  in  direction  ( 9r ,  </>,)  owing  to  light 
incident  onto  the  surface  from  all  directions.  We  then  rewrite  (7.5)  as 

dLr  (6>  4)  =  BRDF ( 0t , ,9r,<j>r)L{^ 0t ,  fa )  cos 0t  dQ, 


and  then  integrate  over  all  incident  directions  to  get  the  total  reflected  radiance  in  direction  (0r, 

m  21]: 

Lr(Pr4r)  =  f  A  ( o, , $ )  BRDF ( 0, , 0 , 9r , £ )  cos tfy/Q, 

2  71 : 

(7 

s  f  L,{0nfi)r(0»0n9r’4r)dnr 


where  r(0„<f>l,9r,</)r)  is  the  radiance  reflectance  function.  Clearly,  r( 9„ fc, 0r, <j>r)  =  cos$  BRDF (Ot,&i,0n&r), 
and  the  two  functions  are  equivalent  ways  of  describing  a  surface.  Some  definitions  of  BRDF  include  a 
factor  of  n  in  the  numerator  of  Eq.  (7.5).  Note  that  the  BRDF  is  a  reflectance  per  unit  solid  angle;  it  can 
have  any  non-negative  value.  As  we  will  see  below,  only  when  the  BRDF  is  integrated  over  solid  angle 
to  get,  for  example,  an  irradiance  reflectance  that  the  result  is  bounded  by  one. 

It  is  emphasized  that  the  BRDF  completely  describes  the  net  effect  of  everything  that  happens  on 
or  below  the  surface  where  it  is  measured.  For  example,  if  the  BRDF  is  measured  in  the  water  column  1 
meter  above  a  sea  grass  bed,  then  all  the  effects  of  the  light  interacting  with  the  grass,  sediments,  and 
water  below  the  1  meter  surface  are  accounted  for  in  this  BRDF.  Knowing  the  BRDF  on  this  imaginary 
surface  would,  for  example,  allow  Hydrolight  (Sequoia  Scientific)  to  compute  the  radiance  distribution  in 
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the  region  above  the  depth  where  the  BRDF  was  measured.  Predicting  or  computing  the  BRDF  of  the 
grass  and  sediments  is  quite  another  story:  to  do  that  one  must  understand  all  of  the  extremely 
complicated  interactions  of  light  with  the  grass  and  sediment  particles. 

In  a  Monte  Carlo  simulation,  we  track  many  individual  photon  packets  as  they  interact  with  the 
medium  and  its  boundary  surfaces.  In  this  case,  the  BRDF  must  be  used  as  a  probability  distribution 
function  to  determine  the  direction  and  weight  of  the  reflected  photon  packet  whenever  a  photon  packet 
hits  the  bottom.  Given  a  photon  packet  with  weight  vv,  that  is  incident  onto  the  bottom  in  direction  ( 0„  f) 
we  need  to  compute  the  weight  wr  and  direction  ( 9r,  (/>,  )  of  the  reflected  photon  packet. 

Since  the  input  direction  ( 9„ (f>,)  is  known,  the  BRDF (0i,</)i,0r,</)l)  can  be  viewed  as  an 
(unnormalized)  bivariate  PDF  for  the  reflected  angles  9r  and  (j)r.  Note  that,  in  general,  these  angles  are 
correlated.  The  directional-hemispherical  reflectance  pdh  for  the  given  (#,$)  can  be  computed  with 


=  J  BRDF  (9^,9^)  cos9r  dil(Orjf)r) 

27T/ 

2k  k  (7-7) 

=  J  J  BRDF  ( 9i  ,t/)l,9r,(j)r)  cos  6r  sin  0r  d  0r  d  <f>r . 

0  0 


The  reflected  photon  weight  is  reduced  by  pdh  (i.e.,  wr  =  pdh( 0„  f)  w, ).  To  determine  the  new  photon 
direction,  we  must  compute  the  cumulative  distribution  functions  for  f  and  for  0,.  The  CDF  for  the 
azimuthal  angle  is  given  by 


-t  <t>r  Till 

CDFMr)=  dh  ( n  ,L  \  1  I  BRDF (Q,  ft, 9, </>)cos9sin9 d9 d</> . 
P  o 


(7.8) 


Note  that  the  directional-hemispherical  reflectance  is  being  used  to  convert  the  BRDF  into  a  normalized 
bivariate  PDF  for  9,-  and  (/>,.  We  are  then  integrating  out  the  0,  dependence  to  leave  a  PDF  for  tj)r ,  which  is 
then  being  used  to  construct  the  CDF  for  (j>r.  Setting  the  CDF  equal  to  a  random  number  p  from  a  uniform 
[0,1]  distribution,  we  must  solve  the  equation 


P  =  CDFt(fr)  (7.9) 

for  tj)r  to  obtain  the  randomly  determined  azimuthal  angle  of  the  reflected  photon  packet.  The  CDF  for  9r 
is  given  by 


CDF,(0r) 


A 

J  BRDF(Oi ,  (j)n0,  (f)r )  cos  #sin  Odd 

_o _ 

n!  2 

J  BRDF (9n 9, (j)r) cos 9 sin 9 d9 

o 


(7.10) 
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Note  that  the  angle  <j>r  determined  in  Eq.  (7.9)  is  used  in  the  BRDF  in  Eq.  (7.10)  when  evaluating  the  0 
integrals.  This  is  how  the  correlation  between  and  <j),  is  accounted  for  in  the  determination  of  the 
reflection  angles.  Finally,  we  draw  another  random  value  p  and  obtain  the  polar  angle  of  the  reflected 
photon  packet  by  solving  the  equation 


p  =  CDFs(er)  (7.11) 

for  0r.  For  all  but  the  simplest  BRDFs,  the  above  equations  must  be  evaluated  numerically  for  each 
photon  packet,  which  can  be  an  enormous  computer  cost  when  millions  of  photon  packets  are  being 
traced. 


7.2.2  Lambertian  Surfaces 

In  many  applications  we  can  ignore  the  directional  dependence  of  the  reflection,  which  is  equivalent 
to  assuming  that  the  surface  is  a  Lambertian  surface.  A  Lambertian  surface  by  definition  reflects  radiance 
equally  into  all  directions.  There  is  a  subtlety  in  this  statement.  For  a  given  incident  lighting,  the  number 
of  photons  reflected  by  each  point  of  a  Lambertian  surface  is  proportional  to  cos#.,  which  is  why 
Lambertian  surfaces  are  sometimes  called  “cosine  reflectors.”  Flowever,  if  we  view  the  surface  with  a 
radiance  detector  having  a  fixed  field  of  view,  the  area  of  surface  that  we  are  viewing  is  proportional  to 
1/cos#.  Thus  the  number  of  photons  going  into  the  detector  is  independent  of  Or,  and  the  reflected 
radiance  is  independent  of  direction.  Its  BRDF  is  simply  equal  to  the  constant  Rfu,  where  Rb  is  called  the 
reflectivity  of  the  surface  or,  in  this  case,  the  bottom  albedo.  The  reflectivity  varies  from  zero  for  a 
completely  absorbing  (“black”)  surface,  to  unity  for  a  completely  reflecting  surface.  There  are  no 
Lambertian  surfaces  in  nature,  but  matte  paper  is  a  good  approximation  except  at  grazing  angles  (0,  and  0r 
near  90  degrees),  where  the  surface  begins  to  look  shiny. 

For  a  Lambertian  seafloor  with  albedo  Rb,  we  can  either  let  the  photon  reflect  if  and  only  if  q  <  Rh  or 
reflect  all  photons  and  multiply  the  weight  by  Rb.  In  either  case,  the  new  photon  direction  is 

pz  =  -yfq. 


(j)  —  2  k  q, 

Mx  =  ^-M2zcos</>, 

My  =  VwJsinf 

(Recall  that  each  appearance  of  q  is  an  independent  random  number.)  The  new  depth  z'  for  a  photon 
reflected  off  the  bottom  is 


z'  =  -  (z-zb), 


where  z  is  the  depth  (z  >  zb)  the  photon  would  have  reached  in  the  absence  of  the  bottom. 
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7.2.3  Minnaert  BRDF 

The  Minnaert  BRDF  is 


BRDF ( 6* , tf), ,0r, </>,.)  =  —( cos 6*  cos 6r . 

Note  that  for  k  =  0  this  reduces  to  the  Lambertian  BRDF.  This  BRDF  was  created  to  explain  the  curious 
fact  that  the  full  moon  appears  almost  uniformly  bright  from  the  center  to  the  edge  of  the  lunar  disk.  If 
the  lunar  dust  were  a  Lambertian  reflector,  the  full  moon  would  appear  bright  at  the  center  and  darker  at 
the  edge.  However,  the  Minnaert  BRDF  agrees  with  observation  over  only  a  limited  range  of  angles. 

Equations  (7.7)  to  (7.11)  can  be  evaluated  analytically  for  the  Minnaert  BRDF.  Equation  (7.7) 
evaluates  to 


dh 

P  = 


2p-cosk  a. 


k  +  2 


which  reduces  to  pAh  =p  for  a  Lambertian  surface.  Equation  (7.8)  gives  just 


CDFt(l) 


2  n 


Plugging  this  into  Eq.  (7.9)  and  solving  for  (j>r  gives 

<j)r  =  2  np  . 

Thus  the  azimuthal  angle  is  uniformly  distributed  over  2 n  radians.  The  CDF  for  0r  as  given  by  Eq.  (7. 10) 
is 


P(6>)  =  l-cosA+16>  . 


Equation  (7.1 1)  then  gives 


Or  =  COS"1  ((/>)' ^+2)) 

after  noting  that  1  —  p  has  the  same  distribution  as  p.  For  a  Lambertian  surface,  the  randomly  generated 
0r  angles  are  distributed  as  cos'1.!/?12),  which  certainly  isn’t  intuitive.  However,  this  distribution  is 
precisely  what  is  necessary  to  make  the  number  of  reflected  photons  per  unit  solid  angle  proportional  to 
cos  0,. 
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8.  BACKWARD  MONTE  CARLO  SIMULATIONS 

When  modeling  three-dimensional  ocean-atmosphere  problems,  it  is  often  necessary  to  use  Backward 
Monte  Carlo  simulations.  The  Backward  Monte  Carlo  (BMC)  approach  is  most  useful  for  an  extended 
source  (e.g.,  sky  radiance  onto  the  sea  surface)  and  a  point  (or  small)  detector,  which  is  what  commonly 
exists  in  oceanography.  Only  BMC  lets  us  simulate  a  point-sized  detector.  Backward  Monte  Carlo 
simulation  generally  refers  to  ray  tracing  of  photons  in  the  reverse  direction;  that  is,  to  trace  photons  from 
the  detector  to  the  source  rather  than  from  the  source  to  the  detector.  Backward  MC  is  useful  because  it  is 
wasteful  to  follow  simulated  photons  in  the  forward  direction  when  very  few  that  are  incident  on  the  sea 
surface  are  actually  collected  by  the  simulated  detector.  Backward  techniques  allow  us  to  only  follow 
photons  that  are  pertinent  to  the  final  outcome  of  the  simulation.  In  addition,  in  the  forward  direction  it  is 
impossible  to  know  how  large  an  area  of  the  sea  surface  should  receive  incident  photons  because  that  area 
depends  on  the  local  factors  of  sky  conditions,  water  lOP’s,  and  the  location  and  orientation  of  the 
simulated  instrument. 

Here  we  will  only  discuss  backward  Monte  Carlo  simulations  in  the  context  of  problems  consisting 
of  a  body  of  water  that  is  illuminated  at  the  sea  surface.  To  implement  a  BMC  simulation,  we  initiate 
photons  at  the  position  of  the  detector.  The  photon’s  initial  direction  is  determined  by  sampling  from  a 
cumulative  distribution  function  (CDF)  designed  to  reflect  the  radiance  response  of  the  detector.  Photons 
that  reach  the  air-side  of  the  air-water  interface  are  weighted  by  the  probability  that  a  photon  would  be 
incident  on  the  seasurface  in  the  exact  opposite  direction  for  the  given  illumination  conditions. 

For  example,  consider  an  upward-facing  irradiance  detector  that  measures  Ed.  Because  an  irradiance 
detector  has  a  cosine  response  to  radiance,  i.e., 


Ed  =  f  M,  L(M)  dju , 

o 

the  probability  density  function  (PDF)  for  emitting  photons  from  an  irradiance  sensor  in  a  BMC 
simulation  is  proportional  to  ju:  and  the  CDF  is  proportional  to  ///’, 

n 

f  0  fO  2 

P(u  )  =  2n\  p(ju)djU  =  2  judju  =  -ju  -q  . 

J  JUS  J  Us 

We  therefore  generate  initial  photon  directions  for  a  BMC  simulation  with 

/T  =  ~yfq  , 


O  -2k q. 

Shown  in  Table  6  are  example  functions  for  generating  the  initial  directions  for  a  variety  of  ideal 
sensors.  As  always  in  this  report,  ju:  =  1  is  taken  to  be  the  downward  direction. 
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Table  6  -  Generating  BMC  functions  for  various  detector  types. 


forward-problem 
detector  type 

backward  MC  source 
function 

Ed 

=  sfq 

Eu 

>T  =s[d 

Eod 

M,  =  -Q 

Eou 

H-.  =  7 

E0 

J: 

ll 

l 

ro 

-o 

cosine  response 
sensor 

^  =V1“9,sin  2df 

For  non-ideal  sensors,  a  CDF  can  be  constructed  using  laboratory  measurements.  For  example, 
for  the  Flyper-TSRB  (Satlantic,  Inc.),  the  CDF  is  given  in  Figure  9  [3,25], 


0  5  10  ,  ,  15  20  25 

angle  (degrees) 

Figure  9  -  Cumulative  distribution  function  for  the  Hyper-TSRB  upwelling  radiometer  (derived  from  data  provided  by  Satlantic). 


See  Refs.  7  and  21  for  more  detailed  discussions  of  the  theory  behind  Backward  Monte  Carlo 
calculations. 


9.  SUMMARY  OF  FREQUENTLY-USED  MONTE  CARLO  EQUATIONS 

Many  practical  ocean  optics  Monte  Carlo  simulations  can  be  formed  with  the  following  short  list  of 
equations. 

The  photon  geometric  pathlength  is  computed  with 

l  ln(<7) 

s  =  —  = - — . 

c  c 

Given  this  pathlength,  the  photon  is  moved  in  the  direction  of  its  direction  cosines  with 
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V 

X 

M>: 

y' 

- 

y 

+ 

My 

z' 

z 

_M2_ 

The  azimuthal  scattering  angle  is 


®  =  2  nq  . 

The  polar  scattering  angle  depends  on  the  scattering  phase  function  being  used.  In  some  cases  the  angle 
can  be  expressed  explicitly,  e.g., 


Ms  ~  cos Tj  =\-2q , 


(isotropic  scattering) 


Ms  = 


2  g 


\  +  g2- 


1-g 


X2 


\  +  g-2gR 


,  g  ^  0  .  (Henyey-Greenstein  scattering) 


In  other  cases,  the  polar  scattering  angle  must  be  interpolated  from  tabulated  values.  The  direction  cosines 
can  be  updated  with 


m'x 

MxMz/^-mI 

1 

1 

J: 

N  tO 

Mx 

\ll~Ms  cos® 

My 

= 

MyM-J^-Ml 

mJ^-mI 

My 

x/l  “  Ms  sin  ® 

A. 

-•fM 

0 

M: 

Ms 

Ms 

yJl~Ms  COS® 

My 

ii 

^1  -//s2  sin  ® 

,  M>  1 

Az. 

Ms 

When  a  photon  reaches  an  air-water  interface  from  the  water  side,  the  incident,  reflected,  and 
transmitted  angles  relative  to  the  surface  are 


9t  -  cos  1  |  ju_  ■  h  | 


6,=0„ 

9,  =  sin-1 


—sin  9. 

\ni  j 


and  the  fraction  of  reflected  light  is  given  by 
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p(0,A)=\ 


sin(t 

%-o,) 

2 

+ 

tan(< 

%-e 

) 

1 

sin(i 

?/  +  3)_ 

tan  (t 

%  +  e 

)J 

l 

p{6i)  = 


\2 


«j  -  n, 

V  77,  +  77zy 


,  0,=o. 


When  a  photon  reaches  an  air-water  interface  from  the  water  side,  the  angles  of  interest  are 

Oj  =  cos-1  /4  •  77  | . 


6,  =  sin 


—sin  0t 

\n<  J 


The  fraction  of  reflected  light  is  given  by 


7>W.0,)  =  | 


sin(t/  -  Qt ) 
sin((9,  +  9t ) 


+ 


tan  (3  -Of) 
tan  ($,  +  6t ) 


l  C7  l 


where 


p{th,n2): 


(  77 j  -  772  ^ 
V  W1  +  "2  y 


,  0,=O 


M3.3)  =  i.  3^3 


3.  =  sin'1 


^  77,  ^ 


Reflection  off  a  Lambertian  surface  is  handled  with 

Ms  =-'IcT  (f)  =  27tq. 


10.  EXAMPLE  APPLICATIONS 

10.1  Collimated  Light  In  A  Non-Scattering  Medium 

Shown  in  Figure  10  is  a  simple  Matlab  script  for  computing  the  light  intensity  at  a  specified  distance 
from  a  collimated  source  in  a  non-scattering  medium.  Photons  are  emitted  at  position  z  =  0,  and  the 
distance  they  travel  (s)  is  determined  with  Eq.  (3.4)  for  the  specified  attenuation  coefficient  (c)  and  the 
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random  number  returned  by  the  Matlab  function  rand.  Photons  that  pass  the  detector  position  (zd)  are 
counted.  The  total  measured  intensity  E  is  expressed  as  a  fraction  of  the  source  intensity.  The  theoretical 
result  is 


E  =  exp(-czd), 

which,  for  c  =  0.5  and  zd  =  3,  yields  E  =  0.2231.  In  a  test  of  the  code,  running  the  script  in  Figure  10  ten 
times  with  N=  105  photons  returned  a  mean  value  of  0.2230  with  a  standard  deviation  of  0.0014. 


N  =  le5;  %  number  of  photons 
zd  =  3;  %  position  of  detector 

c  =  0.5;  %  attenuation  coefficient 

E  =  0;  %  initialize  detector  response 

for  i=l:N 

z  =  0;  %  initial  photon  position 

s  =  -log(rand)/c;  %  geometric  path  length 
z  =  z  +  s;  %  move  photon 

if  z  >  zd,  E=E+1 ;  end  %  count  photons  that  reach  detector 
end 

E  =  E/N _ %  normalize  detector  response _ 

Figure  10  -  Matlab  Monte  Carlo  code  for  simulating  collimated  light  in  a  non-scattering  medium. 

10.2  Irradiances  Out  The  Top  Of  Scattering  Medium 

Shown  in  Figure  11  is  a  simple  Monte  Carlo  Matlab  script  for  computing  the  intensities  returned 
from  an  isotropically  scattering  semi-infinite  medium.  The  slab  is  illuminated  with  collimated  light  at 
polar  angle  cos_1(//z).  The  distance  each  photon  travels  before  its  first  interaction  is  determined  by  Eq. 
(3.4),  at  which  point  the  photons’s  weight  is  multiplied  by  the  single-scattering  albedo  oy  and  the  photon 
is  scattered  isotropically  [Eq.  (4.11)].  If  and  when  a  photon  exits  through  the  top  of  the  scattering 
medium,  then  it’s  weight  is  added  to  the  count.  The  count  normalized  by  the  number  of  photons  traced 
gives  the  intensities  leaving  the  top  of  the  medium.  Benchmark  numerical  results  for  this  problem  are 
provided  in  Van  de  Elulst  [31,  Table  12].  Example  results  from  our  Matlab  script,  along  with  the 
corresponding  Van  de  Flulst  results,  are  given  in  Table  7. 
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c  =  1;  %  attenuation  coefficient  (1/m) 

wO  =  0.8;  %  single-scattering  albedo 

N  =  le5;  %  number  of  photons  to  trace 

ns  =  100;  %  max  number  of  scatters  per  photon 

E  =  0;  %  initialize  detector 

for  i=l:N 


z  =  0;  %  initial  position  depth 

muz  =  0.1;  %  incident  light  direction  (collimated) 

w  =  1 ;  %  initial  photon  weight 

for  j=l:ns 

s  =  -  log(rand)/c;  %  geometric  path  length 

z  =  z  +  muz*s;  %  move  photon 

if  (z<0)  E=E+w;  break,  end  %  count  photons  leaving  out  top 
w=w*w0;  %  absorb  fraction  of  photon  packet 

muz=  1  -  2* rand;  %  isotopic  scattering 

end 
end 

E  =  E/N  %  normalize  result 


Figure  1 1  -  Monte  Carlo  Matlab  script  for  computing  the  intensities  returned  from  an  isotoprically  scattering  semi- infinite 

medium. 


Table  7  -  Comparisons  between  results  from  code  in  Figure  1 1  with  those  from  Van  de  Flulst  [31]. 


A) 

Wo 

Ref.  31  Eu(0) 

MC  Eu( 0)  mean  ;  standard  deviation 
( 10  runs  with  jV=105  each) 

1.0 

0.8 

0.28525 

0.2854255  ;  0.000811 

1.0 

0.4 

0.08336 

0.0835981  ;  0.000350 

0.1 

0.8 

0.49071 

0.4907184;  0.001170 

10.3  Updating  Direction  Cosines 

Shown  in  Figure  12  is  a  Matlab  function  for  updating  the  direction  cosines  during  a  scattering  event. 
This  is  necessary  for  simulating  anisotropic  scattering,  such  as  that  in  water. 
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function  [alp, bet, gam]  =  chgdir( alpha,  beta,  gamma,  gammas,  phis) 

%  Determines  the  new  direction  cosines  (alp, bet, gam)  from  old 
%  direction  cosines  (alpha,  beta,  gamma),  polar  scattering  cosine 
%  (gammas),  and  scattering  azimuthal  angle  (phis). 

% 

%  Robert  A  Leathers  &  Trijntje  Downes,  17  March  2000. 

stheta  =  sqrt(l-gamma*gamma);  %  sin(theta)  of  previous  photon  direction 

sthetas  =  sqrt(l-gammas*gammas);  %  sin(theta)  of  scattering  angle 
alphas  =  sthetas*cos(phis);  %  x-direction  cosine  of  scattering 

betas  =  sthetas*sin(phis);  %  y— direction  cosine  of  scattering 

%  new  photon  directions... 

if  (stheta  >  le-12)  %  if  initial  direction  not  straight  up  or  straight  down 
B  =  [alpha*gamma/stheta,  -beta/stheta,  alpha;... 
beta*gamma/stheta,  alpha/stheta,  beta;... 

-stheta,  0,  gamma]; 

[alphas;betas;gammas]; 

newdir  =  B*[alphas;betas;gammas]; 

alp  =  newdir(l); 

bet  =  newdir(2); 

gam  =  newdir(3); 

else  %  initial  direction  straight  up  or  down  (sin(theta)=0) 

%  then  new  direction  cosines  =  +/-  scattering  direction  cosines 
s  =  sign(gamma); 
alp  =  s*  alphas; 
gam  =  s*gammas; 
bet  =  s*betas; 

End _ 

Figure  12  -  Matlab  function  for  updating  the  direction  cosines  for  anisotropic  scattering. 

The  function  in  Figure  12  was  used  to  generate  Figure  13.  An  initial  photon  direction  was  selected  to 
be  jux  =  0.2,  juy  =  -  0.1,  and  /jx  =  sqrt(l  -  /A  -  /A).  The  scattering  polar  angle  was  selected  such  that 
cosxF  =  0.4,  and  the  scattering  azimuthal  angle  was  selected  randomly  using  Eq.  (4.3).  Shown  in  Figure 
13  are  the  endpoints  of  100  unit  vectors  computed  with  function  chgdir.  It  can  be  seen  that  Eq.  (4.3) 
provides  a  uniform  distribution  of  the  polar  scattering  angle. 
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10.4  3-D  Diffusion  From  Point  Source 

The  code  in  Figure  14  computes  the  absorption  at  a  specified  radius  away  from  a  point  source  in  an 
isotropically  scattering  medium. 


a  =  1 ;  %  absorption  coefficient  (1/m) 
c  =  6;  %  attenuation  coefficient  ( 1/m) 
wO  =  (c-a)/c;  %  single-scattering  albedo 
P0  =  1 ;  %  source  intensity  (W) 

N  =  1  e4;  %  number  of  photons  to  trace 

ns  =  30;  %  number  of  scatters  to  trace 
Nd  =  20;  %  number  of  detectors 

dd  =  0.02;  %  detector  spacing  (m) 

A  =  zeros(l,Nd);  %  initalize  detectors  (all  photons) 
for  i=l:N 

p=[0,0,0];  %  initial  position  (x,y,z) 

muz  =  l-2*rand;  %  isotropic  source 
phi  =  2*pi*rand;  %  isotropic  source 

mux  =  sqrt(l-muzA2)*cos(phi);  %  initial  x-direction  cosine 
muy  =  sqrt(l-muzA2)*sin(phi);  %  intitial  y-direction  cosine 
r=0;  %  intial  distance  from  source 
w  =  P0;  %  initial  photon  weight 

for  j=l:ns 

s  =  -log(rand)/c;  %  geometric  path  length 

p  =  p  +  s*[mux,muy,muz];  %  move  photon 

r=  norm(p);  %  radial  position _ 
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index  =  ceil(r/dd);  %  find  detector  number 
if  (index<=Nd),  %  is  photon  within  detector  grid? 

A(index)  =  A(index)  +  (l-wO)*w;  %  count  absoiption 
end 

w  =  w*wO; 

mu=  1  -  2* rand;  %  isotopic  scattering 
phis  =  2*pi*rand;  %  isotopic  scattering 

[mux,muy,muz]=chgdir(mux,muy,muz,mu,phis);  %  new  direction 
end 
end 

r  =  [1  :Nd]*dd;  %  radii  for  detectors 
A  =  A/N;  %  normalize  detector  response 
A1=A1/N;  %  normalize  detector  response 
A  =  A./(dd*4*pi*r.A2);  %  convert  response  to  W/mA3 
ht  =  3*a*c*P0*exp(-r*sqrt(3*a*c))./(4*pi*r);  %  theory 

semilogy(r,ht,r,A,'*') 

grid  on 

xlabel('r'),  ylabel('W/mA3') 

legend('theory','MC  (total)', 'MC  (positive  xyz)') _ 

Figure  14  -  Matlab  Monte  Carlo  code  for  computing  the  absorption  as  a  function  of  radius  away  from  a  point  source  in  an 

isotropically  scattering  medium. 


The  theoretical  result  is  given  by 

A  =  ^ac^°  exp(-rV3ac)  (10.1) 

4  nr 

The  results  from  the  MC  code  are  compared  with  the  theoretical  result  in  Figure  15. 


Figure  15  -  Comparison  of  MC  result  from  Figure  14  with  theoretical  result  using  104  photons. 
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10.5  Optically  Deep  Ocean 

Kirk  [9]  provides  example  FORTRAN  code  for  a  simple  1 -dimensional  ocean  model.  Collimated 
light  is  incident  on  a  level  sea  surface.  Refraction  of  the  incident  light  is  included  but  reflection  is  not. 
Internal  reflection  is  treated  as  a  step  function,  where  only  light  below  the  critical  angle  is  reflected. 
Water  scattering  is  described  by  Petzold's  measurements  in  San  Diego  Harbor.  The  seafloor  is  taken  to  be 
100%  absorbing  and  is  intended  to  be  located  optically  far  from  the  surface.  No  photon  weighting  is  used. 
Checks  are  put  in  place  to  measure  radiance  in  5  degree  angular  bins. 

Alternatively,  Figure  16  shows  our  Matlab  code  for  computing  upwelling  and  downwelling 
irradiances  in  a  homogeneous  ocean.  Scattering  is  described  by  the  Henyey-Greestein  scattering  phase 
function.  Results  from  the  ocean  simulation  of  Figure  16  agreed  well  with  those  from  Hydrolight 
(Sequoia  Scientific,  Redmond,  WA). 


c  =  1;  %  attenuation  coefficient  (1/m) 

wO  =  0.8;  %  single-scattering  albedo 

g  =  0.9;  %  scattering  asymmetry  factor 

zd  =  1.0;  %  depth  of  downwelling  detector 

zu  =  1 .0;  %  depth  of  upwelling  detector 

n  =  1 .34;  %  relative  index  of  refraction  for  water 

N  =  le4;  %  number  of  photons  to  trace 

ns  =  50;  %  number  of  scatters  to  trace 

muc  =  sqrt(l-l/nA2);  %  critical  angle  for  internal  reflection 

Eu=0;  Ed=0;  %  initiate  detectors 

for  i=l:N 

z  =  0;  %  initial  position  depth 

muz  =  sqrt(rand);  %  diffuse  illumination 
w  =  1 ;  %  initial  photon  weight 
if  (muz==l), 

ra=  ((n-l)/(n+l))A2;  %  reflection  normal  to  surface 
else 

ti  =  acos(muz);  %  angle  of  incidence 
tt  =  asin(sqrt(l-muzA2)/n);  %  air-to-water  trans.  angle 
ra  =  ((sin(ti-tt)/sin(ti+tt))A2+(tan(ti-tt)/tan(ti+tt))A2)/2; 
muz  =  cos(tt);  %  new  direction 
end 

w  =  w*(l-  ra);  %  air-to-water  transmission 
for  j=l:ns 

zold=z;  %  previous  position 
s  =  -log(rand)/c;  %  geometric  path  length 
z  =  z  +  muz*s;  %  move  photon 
if  (z<zu)&(zold>zu),  %  photon  passes  Eu  sensor 
Eu=Eu+w;  %  update  upwelling  sensor 
end 

if  (z<0)  %  photon  reaches  surface 

muz  =  -  muz;  %  change  to  positive  value 
z  =  -  z;  %  new  position  of  reflected  photon 
zold  =  0;  %  photon  now  traveling  down 
_ if  (muz>muc),  %  partial  transmission _ 
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if(muz==l),  %  photon  vertical 

ra=  ((n-l)/(n+l))A2;  %  reflection  normal  to  surface 
else 

ti  =  acos(muz);  %  angle  of  incidence 
tt  =  asin(n*sqrt(l-muzA2));  %  water-to-air  trans.  angle 
rw  =  ((sin(ti-tt)/sin(ti+tt))A2+(tan(ti-tt)/tan(ti+tt))A2)/2; 
w  =  w*rw;  %  weight  of  reflected  photon 
end 
end 
end 

if  (z>zd)&(zold<zd),  %  photon  passes  Ed  sensor 
Ed=Ed+w;  %  update  downwelling  sensor 
end 

w=w*wO;  %  in-water  scattering 

mu  =  (l+gA2-((l-gA2)/(l+g-2*g*rand))A2)/(2*g);  %  HG  scattering 
Phi  =  2*pi*rand;  %  azimuthal  scattering  angle 

muz  =  muz*mu  -  sqrt(l-muzA2)*sqrt(l-muA2)*cos(Phi);  %  new  direction 
end 
end 

disp(sprintf('Eu  at  %g  m:  %g',zu,Eu/N)); 

disp(sprintf('Ed  at  %g  m:  %g',zd,Ed/N)); _ 

Figure  16  -  Matlab  code  for  the  MC  simulation  of  a  homogeneous  water  column. 


10.6  PSICAM  Analysis 

In  the  analysis  of  a  Point-Source  Integrating-Cavity  Absoiption  Meter  (PSICAM)  [1],  it  is  useful 
to  know  the  probability  that  a  photon  leaving  the  inner  wall  of  cavity  will  return  to  the  wall  (as  opposed  to 
being  absorbed  by  the  water  within  the  cavity).  For  given  water  optical  properties,  this  value  can  be 
computed  with  the  Matlab  code  provided  in  Figure  17.  In  the  absence  of  scattering,  the  theoretical  result 
is 


P  = 


l-(2ar  +  l)exp(-2ar) 


2a2 


where  r  is  the  inner  radius  of  the  cavity.  The  provided  code  can  be  used  to  compute  the  value  for  any 
amount  of  scattering.  It  can  also  be  modified  to  account  for  non-ideal  characteristics  of  a  particular 
instrument. 
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a  =  0.2;  %absorption  coeff 

b  =  50;  %scattering  coeff 

model  =  1;  %  l=isotropic,  2=HG,  3=Petzold 

g  =  0.92;  %scattering  asymmetry  factor  (required  for  HG  only) 

r  =  0.05;  %  radius  of  sphere  (m) 

N=1000;  %number  of  photons  traced 

minweight=0.0001;  %  min  photon  weight  before  discarding  it 
maxns=10;  %maximum  number  of  scatters  for  each  photon 

c  =  a+b;  %  attenuation  coefficient 
wO  =  b/c;  %single-scattering  albedo 

count=0; 

for  i=l  :N,  %  trace  photons  one  at  a  time 

weight  =  1 ;  %  new  photon,  generated  at  position  (0,0, -r):  bottom  of  sphere 

ns  =  0; 
y  =  0; 
z  =  -r; 
x  =  0; 

val=(i-0.5)/N;  %step  uniformly  over  interval  (0,1) 
mu=sqrt(1.0-val);  %  diffuse  light  directions 

gamma  =  mu;  %  alpha,  beta,  gamma  =  direction  cosines 
beta  =  0;  %move  in  x-z  plane 

alpha  =  sqrt(1.0-muA2);  %alphaA2+betaA2+gammaA2=l 

while  (ns  <=  maxns)  &  (weight  >  minweight)  %  keep  tracing  photon  path 
%how  far  does  it  go? 
l=-(l/c)*log(l-rand);  %totalpath  length 

z  =  z+gamma*l;  %new  position 
x  =  x+alpha*l; 
y  =  y+beta*l; 

if  ((xA2+yA2+zA2)  >  rA2);  %does  it  reach  the  wall? 
count=count  +  weight; 

break  %if  the  photon  hits  the  wall,  break  out  of  while  loop  and  don't  scatter  it. 
end 

%  scatter  photons 

weight=weight*w0;  %  fraction  scattered  (versus  absorbed)  given  by  wO 
ns  =  ns+1;  %number  of  times  the  photon  has  been  scattered 
gammas  =  scat(model,g);  %  compute  polar  scattering  angle 
phis  =  2*pi*(rand); 

[alpha.beta, gamma]  =  chgdir(alpha, beta, gamma, gammas, phis); 
end 
end 

F=count/N; 

disp(sprintf('The  probability  of  photon  survival  from  wall  to  wall  is  %g',F)) 


Figure  17  -  Matlab  code  for  computing  the  probability  of  photon  survival  from  wall  to  wall  within  a  PSICAM. 
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